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Stochastic aspects of chemical reaction models related to the Soai reactions as well as to 
the homochirality in life are studied analytically and numerically by the use of the master 
equation and random walk model. For systems with a recycling process, a unique final prob- 
ability distribution is obtained by means of detailed balance conditions. With a nonlinear 
autocatalysis the distribution has a double-peak structure, indicating the chiral symmetry 
breaking. This problem is further analyzed by examining eigenvalues and cigenfunctions of 
the master equation. In the case without recycling process, final probability distributions 
depend on the initial conditions. In the nonlinear autocatalytic case, time-evolution starting 
from a complete achiral state leads to a final distribution which differs from that deduced 
from the nonzero recycling result. This is due to the absence of the detailed balance, and a 
directed random walk model is shown to give the correct final profile. When the nonlinear 
autocatalysis is sufficiently strong and the initial state is achiral, the final probability dis- 
tribution has a double-peak structure, related to the enantiomeric excess amplification. It 
is argued that with autocatalyses and a very small but nonzero spontaneous production, a 
single mother scenario could be a main mechanism to produce the homochirality 
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1. Introduction 

For some organic molecules two kind of stereostructures that are mutually mirror sym- 
metric are possible to exist, and they are called enantiomers. 1 Two enantiomers are chiral 
like the left- and the right-hand which cannot be overlapped by translational and rotational 
transformations. Since their physical properties are mostly the same except the response to 
the optical polarity, there should be an equal amount of both enantiomers in the production 
started from achiral substrates. In nature, on the other hand, it has long been known that the 
chiral symmetry in life is broken; 2,3 all proteins consist of the left-handed L-amino acids and 
nearly all sugars belong to the right-handed D-series. 4 The origin of this homochirality has 
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attracted much attention in relation to the origin of life itself 1, 3,5-7 Ideas explaining the origin 
of homochirality are categorized in two groups; it is brought by some external advantage factor 
for one chirality to the other, 3 or it happens by accident. 8 In both cases, however, the degree 
of excess of one enantiomer to the other is expected very small, 6,9 and the amplification of 
enantiomeric excess (ee) is indispensable. 

In an open system, Frank proposed long ago a theoretical model which contains auto- 
catalysis and an antagonistic nonlinear process, and he showed that the model leads to a 
chirality selection. 10 Recently Soai and his coworkers found experimental systems which show 
the ee amplification. 11,12 Experiments were performed in a closed system, and the ee amplifi- 
cation was shown to depend on the initial condition. This result is explained by assuming the 
quadratic autocatalysis. 13,14 By including a recycling process in addition, it is shown that the 
system relaxes to a unique final state with a broken chiral symmetry. 15 

Most of the theoretical analyses have been performed in terms of deterministic rate equa- 
tions. 13-24 In the meanwhile, sequences of experimental runs starting from an initial state with 
no chiral ingredients have been performed. 25-28 In some runs, the system has a preference to 
one chirality, and in some other runs to the other chirality. Values of an ee order parameter are 
distributed wide in a whole possible range. The probability distribution for the ee value has 
symmetric double peaks at intermediate values of ee with opposite signs. Symmetric profile 
reflects the fact that there is no preference in the chirality in the initial situation, and the 
double-peak structure represents that the symmetry breaking is induced in the dynamic evo- 
lution. For the description of the probability distribution of populations of chiral species, we 
have to study the stochastic aspects of the chiral molecule production. There is a stochastic 
study of a spontaneous and a linear autocatalytic chemical reaction without recycling. 29,30 
Here we study stochastic evolution of systems with more generic chemical reactions; linear as 
well as nonlinear autocatalysis with and without recycling. 

In §2, we present a stochastic model for the production of chiral molecules R and S from 
an achiral substrate A in a closed system. Transition probabilities contain processes as a spon- 
taneous production, a linear and a quadratic autocatalytic production, and a back reaction 
which recycles the substrate from the products. In §3, we discuss the shape of the final proba- 
bility distribution by assuming a detailed balance between production and recycling reactions. 
This analysis is valid if the final probability distribution is unique, independent of the inter- 
mediate path. The peak position of the final probability distribution is found to correspond to 
the fixed points of the rate equations. With a quadratic autocatalysis, the distribution has a 
double-peak profile, indicating the chiral symmetry breaking. In §4 the master equation which 
governs the time evolution of the probability distribution is integrated numerically, and the 
final distribution is compared with that obtained in the previous section. In most cases two 
results agree, but if the system has a quadratic autocatalysis without recycling back reaction, 
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the numerical integration gives rise to a profile different from that expected in the limit of 
weak recycling. The chiral symmetry breaking in the case with recycling is interpreted in terms 
of a degeneracy in the eigenvalues of the master equation evolution in §5. The zero eigenvalue 
state corresponds to the final probability distribution, and with a quadratic autocatalysis the 
first non-zero eigenvalue approaches zero in the large number limit of the reactants. With- 
out recycling, there are infinitely many degeneracy of eigenvalues at zero, representing the 
non-ergodicity such that the final probability depends on the initial condition. In §6 a toy 
model for the system without recycling, a directed random walk model, is proposed. It gives 
the final probability distribution for a spontaneous and a linear autocatalytic system analyti- 
cally, and for a quadratic autocatalytic system numerically. The resulting distributions agree 
with those given by the numerical integration of the master equation in §4. The stochastic 
approach reveals (i) the single mother scenario of the homochirality: In the case of a very 
small spontaneous production rate, the first mother chiral species produced by this sponta- 
neous production converts all the substrate molecules into her type of chiral products by some 
autocatalytic processes, before the second mother of another chiral type is born, (ii) With a 
quadratic autocatalysis with a moderate value of the rate constant, double peaks appear in 
a probability distribution at a finite values of the ee order parameter, in agreement with the 
experimental observation: The spontaneously produced racemic single peak splits due to the 
nonlinear amplification of chiral imbalance. The result is summarized in the last section. 

2. Model and Elementary Processes 

In order to understand the stochastic features of chiral symmetry breaking, we use here 
the simplest model considered previously: 15 An achiral substrate molecule A turns into one 
of the two enantiomers R or S in a closed system. Only with spontaneous reactions 

A^R, A^S (1) 

with the same reaction rate ko, it is easily shown that an initial enantiomeric excess (ee) 
decreases. 15 For the ee amplification, some autocatalytic processes are necessary. We con- 
sider two types of autocatalyses; a linear and a quadratic ones. A rate constant of linear 
autocatalyses 

A + R^2R, A + S^2S (2) 

is set ki , and that of quadratic autocatalyses 

A + 2R^3R, A + 2S -» 3S (3) 

is set &2- It is found 15 that additional back reactions from the chiral products R or S to the 
achiral substrate A, which we call recycling processes hereafter, 

R -> A, S -> A (4) 
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select a unique final state with a definite value of the ee: The situation is similar to the sym- 
metry breaking in equilibrium statistical physics, and we call it a chiral symmetry breaking. 
By denoting the rate of the recycling (4) as A, the rate equations of the concentrations of 
chiral enantiomers r and s are written as 

-t: = (&0 + k\r + k 2 r 2 )a - Xr, 

= (ko + kis + k 2 s 2 )a - As. (5) 

Because the system is assumed to be closed, the total concentration c of all the reactant species 
is conserved and the concentration of achiral substrate a is determined as a = c — r — s. The 
order parameter of ee is defined as usual as 

r — s 



(6) 



r + s 

and its absolute value \(j>\ is often called an ee parameter. 

Rate equations in general describe averaged behaviors of the reaction system when there 
are ample amount of molecules. In the initial stage of reaction, however, there are only a very 
few product molecules R and S and the discrete and the stochastic aspect of the chemical 
reaction could be important for autocatalytic reactions. The state of a system is described by 
population numbers of achiral and chiral molecules as Na, Nr, Ng where the total number 
of molecules N are fixed constant so that Na = N — Nr — Ng. At a time t, the system is in a 
state (Na,Nr, Ng) with a probability P(Na, Nr, Ns',t), and the chemical reaction changes 
the state stochastically. 

Each molecule reacts to change its state with a certain transition probability. Let us 
consider first a process of R production such that the state (Na, Nr, Ng) changes to (Na — 
1, Nr + 1, Ng). Since the change is induced by the spontaneous reaction of one A molecule 
among Na of them to R, the transition probability of state per time is given by ^Na- If the 
A molecule encounters one R, the linear autocatalysis increases the reaction rate by k\. When 
there are Nr number of R molecules in a well homogenized volume V, the probability of the 
encounter is equal to the concentration r = Nr/V . Therefore, the increment of the transition 
probability is ki(Nn/V) x Na- With a similar consideration, the quadratic autocatalysis 
increases the transition probability by ^(Nr/V) 2 x A^4. Thus, the total transition probability 
of creating one R molecule is denoted by 

W(N A , N R , N s ^Na- 1, N R + 1, Ng) = (k + kiN R /V + k 2 N 2 R /V 2 )N A 

= (k + k\Nr + k 2 N 2 r )N a (7) 

where we define stochastic rate coefficients 29 ' 30 as 

Kx = ki/V = hc/N, k 2 = k 2 /V 2 = k 2 c 2 /N 2 (8) 
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with c = N/V being the total concentration of chemically relevant molecules. For the recycling 
process (4) with the rate constant A, the transition probability is written as 

W(N A , N R , N S ^N A + 1, N R - 1, N S ) = XN R . (9) 

For the S enantiomer, similar transition probabilities are defined. The time evolution of the 
probability distribution is then written by the master equation as 



dP(N A ,N R ,N s ;t) = £ p^^^w^^ 



s ^N a ,Nr,N s ) 



- P(N A ,N R ,N s ;t)W(N A ,N R ,Ns ^ N' A ,N' R ,N' S ) (10) 

N' R ,N' S 

where the state (N' A , N' R , N' s ) is connected to the state (N A , Nr, Ng) by the transition 
probability W, and differs the latter by numbers of molecules at most by one. Average con- 
centration of R enantiomer is defined by using the probability distribution as 

(r(t)) = (N R (t))/V = V- 1 £ N R P(N A ,N R ,N s ;t) (11) 

N R ,N S 

and its dynamics is described as 

= ((k + hr + k 2 r 2 )a) - X(r(t)). (12) 

Therefore, by neglecting the correlation as (ra) = (r)(a) and (r 2 a) = (r) 2 (a), one recovers the 
usual rate equations (5). 

3. Final probability 

We know that with the recycling process (4), the rate equations (5) lead to a unique 
final state. 15 Therefore, one may expect the existence of a definite probability distribution Pf 
associated with this unique final state. This Pf can be obtained by numerically simulating the 
time evolution (10), but one can obtain it analytically by assuming a detailed balance condition 
such that the production flow from the state (N A + 1, Nr — 1, Ns) to (N A , Nr, N$) written 
as P f (N A + 1,N R - 1, N S )W(N A + l,N R - 1, N s -> N A ,N R , N s ) balances with the counter 
recycling flow P f (N A , N R , N S )W(N A , Nr, N s -» N A + 1, N R - 1, N s ) for N R > 1. This leads 
to the final probability as 

P(N N AT ^ W(N a + 1,Nr-1,N s ^N a ,Nr,N s ) d ^ T AT 

Pf[N A , Nr, Ns) = TT7- / T\r ^7 ITr ITr , -, Ar , A r \ Pf [N A + I, N R - l,N S ) 

W(N A , N R , N s -> N A + 1, N R - 1,N S ) 

A1\R 

N R -1 



(N A + Nr) 



[ (k Q + Kim + K 2 m 2 ) 



P f (N A + N R ,0,N s ). (13) 



N R \N A \ \ N R 
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When N s = 0, Eq.(13) determines Pf(N A ,N R ,0) in terms of P/(iV,0,0). When N s > 1, one 
applies similar procedures for Ns again and obtains the relation 

N R -1 n s -i 

[ (k + Kim + K 2 m 2 ) [ (k + Kin + K 2 n 2 ) 

Pf(N A ,N R ,N s ) = ^U^— xW W<>.°)- 

(14) 

Here the number conservation N A + Nr + Ns = N is used. Pf(N A ,0,Ng) can be obtained 
from Pf(N A ,N R ,0) by simply replacing N R by Ns- The final state probability distribution is 
in fact symmetric for any N R and Ns as 

P f (N A ,N R ,N s ) = P f (N A ,N s ,N R ) (15) 

reflecting the chiral symmetry of the dynamic processes (l)-(4). Pf(N, 0, 0) is determined from 

the normalization condition 

N N-N R 

Yl Pf(N A ,N R ,N s )= Pf(N -N R -N S ,N R ,N S ) = 1. (16) 
N R =0 n s =o n r ,n s 

All the reactions except the recycling (4) consume achiral molecules A. Thus without recycling, 

reactions proceed and come to a halt when all A molecules are consumed, namely A^4 = 0. 

The final distribution Pf should hence be zero for all states with N A ^ as A — > 0. It is 

possible by assuming P/(N, 0, 0) oc A^. Thus the final probability is written as 

P f (N A , N R , N s ) = M NRl N slNA] xNA f( k °' N R )f(ko, ki, k 2 ; N s ) (17) 

with M being the normalization constant, and a function / is defined as 

f(k ,K 1 ,K 2 ;M) = 



1 for M = 0, 

M-l 

J (ko + K\m + K 2 m 2 ) for M > 1. 



(18) 



V m=0 

Its dependence on the stochastic rate coefficients ko, K\ and k 2 is explicitly indicated for later 
convenience. 

For a large system, the probability is denoted as Pf = e~ V9 ^ a ' r ' s ^ with r = N R /V, s = 
Ns/V and a = N A /V = c — r — s, and the effective "potential" g is written as 

g(a, r, s) = - V~ l In P f (N A , N R , N s ) 

= — const + r In r + s In s + a In a — a In A 

— / dx ln(&o + k\x + k 2 x 2 ) — j dyln(ko + k%y + k 2 y 2 ). (19) 
Jo Jo 

The stationary conditions dg/dr = dg/ds = yield 

(ko + k%r + k 2 r 2 )a = Xr, (ko + k\s + k 2 s 2 )a = As, (20) 
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which are identical with the equations to determine fixed points of the rate equations (5). 
There are racemic solutions r = s which satisfy 

(k Q + k 1 r + k 2 r 2 )(c-2r) = Xr (21) 

and chiral solutions r ^ s for positive k 2 as 



with 



'f(W^- 4 ^)= x± - s (22> 

k 2 c 2 - he ± v / {k 2 c 2 + k lC ) 2 - 4k 2 c 2 X 
2k 2 c z 

Chiral state is impossible to exist unless quadratic autocatalysis is active. It is possible only 

when the obtained values of r and s are real and positive. 

The probability around the fixed point is characterized by the second derivatives; 

d 2 g(r, s) 1 1 k\ + 2k 2 r d 2 g(r, s) 1 

Q r 2 r a ^ Q _|_ _|_ ^ 2? .2 ' 9 rs drds a 

_ d 2 g(r,s) = 1 + 1 _ fci + 2k 2 s 
ds 2 s a ko + k\s + k 2 s 2 
Especially around the racemic fixed point, which we denote here r* = s*, symmetry implies 
9rr = 9ss, and the distribution of fluctuations 5r = r — r* and 5s = s — s* is approximately 
given as 

P f {r,s) ^e~ V9 ^ 

* ' 1 

Therefore, the fluctuations are determined with a* = c — r* — s* as 

2 2a*r*(k + k 1 r* + k 2 r* 2 ) 



exp < — V 



-(g rr + g rs )(5r + 5s) 2 + -(g 

rr 9rs 

)(5r-5s) 2 ]y (25) 



{(5r + 5s) 2 ) 



V(g rr + g rs ) V(k c + (c - a*)hr* + (c - 2a*)k 2 r* 2 ) ' 



These formulae are useful in the following detailed studies of the final probability distribution 
for some typical cases. 

3.1 Spontaneous production with and without recycling 

Only with a spontaneous production and a recycling (ko, A > and k\ = k 2 = 0), the final 
probability distribution is easily calculated as 

N R \N S \N A \ (2k + \y 
with N = Na + Nji + Ns- It has a peak at a racemic point obtained from Eq.(21) as 

r* = s*= k ° - c, a * = c-r*-s*= , A - c (28) 



Pf(N A ,N R ,N s ) = ^^^-T^W- ( 2? ) 
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and the probability has a Gaussian profile in the vicinity of the peak as 

P f (N A , Nr, N S ) * exp { - v f^° + A) V + 5s) 2 + ^r±V - 5s) 2 } } (29) 

with 5r = r — r* and 5s = s — s*. The peak position (r*, s*) agrees with the fixed point of the 
rate equations (5). Fluctuations 5r and 5s are anisotropic such that in the direction r = — s 
fluctuation vanishes as A — > whereas in the direction r = s finite fluctuation remains. 

In the limit of vanishing recycling (A — > 0), there is no reason that the final probability 
distribution takes the form given by the detailed balance condition, but it is interesting to 
examine Pf in this limit. Since all the substrate molecules turn into chiral products, Na should 
be zero so that Pf(N^ > 0, Nr, N$) = 0. Thus, the final state lies on a fixed line r + s = c 
or Nr + Ns = N. The probability (27) in this limit takes a binomial distribution along the 
fixed line as 

7V! 

P f (0,N R ,N - Nr) = ; : - . (30) 

A ' K ' HJ N R \(N - Nr)\\2J v ; 

This result is identical with the solution obtained by Lent, who analyzed the same system 

with no recyling assuming a complete achiral state (N, 0, 0) as the initial condition. 29 ' 30 

3.2 Linear autocatalysis with and without recycling 

Addition of a linear autocatalytic process does not alter the final probability distribution 

qualitatively. With k<i = 0, the probability has a peak at a racemic point obtained from 

Eq.(21) as 

*_ V(2fc - he + A) 2 + 8k k lC ' - (2jp - he + A) 
r — s — c — W-U 

Ahc 

which corresponds to a fixed point of the rate equations (5). Fluctuations around the fixed 
point given by Eq.(26) are non-negative as long as &2 = 0. 

If there is no recycling (A = 0), a fixed point is at r* = s* = c/2 and a* = 0. Thus the 
fluctuation ((5r + 5s) 2 } vanishes and r and s fluctuate along the diagonal line r + s = c. 
In this case, fluctuations ((5r) 2 ) = ((5s) 2 ) are still finite as long as ko > 0. Only when the 
spontaneous production is absent (k^ = 0), fluctuations diverge, indicating that a distribution 
becomes flat along a ridge r + s = c. 

In fact, for A — > 0, Eq.(17) tells us that the final probability is nonzero only along the fixed 

line Na = 0, or Nr + iV^ = N, in the Nr-N$ phase space, and the normalized probability is 

calculated explicitly as 

P fn N N at \ AH f(k , Kl ,0;NR)f(k ,K l ,0;N-N R ) 
P f (0,N R ,N-N R ) = NrKn _ Nr)1 f(2k , Kl ,0;N) " (32) 

This result is again identical with the solution obtained by Lent, starting from the complete 

achiral state (N, 0, 0). 29 ' 30 

If ko is as small as ko = k\ = k\/V, then f(ko,ko,0;M) = k^ 1 Ml and f(2ko,ko,0;N) = 

ky (A^+l)!, and the final probability distribution is constant as Pf (0, Nr, N—Nr) = 1/(^+1). 
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For further smaller ko, Pf has double peaks of height Pf = 1/2 at Nr = or N, and 
essentially vanishes otherwise. The peaks at Nr = or N$ = mean that an enantiomer 
first produced by the spontaneous reaction converts all the achiral substrate molecules into 
her type of chirality by linear autocatalysis before the second type of enantiomer is produced 
spontaneously; this is just a single mother scenario. The former process takes place within 
a time about ^jv^=i 1/ k iNr{N — Nr) ~ 2 In N/k\N whereas the latter requires the time 
l/k()N. Therefore, by neglecting the logarithmic correction, the second mother has no chance 
to be born when the spontaneous production rate ko is smaller than n±, 
3.3 Quadratic autocatalysis with and without recycling 

When k 2 > 0, there are at most three racemic fixed points. Those points with r* > s/kojka 
are unstable, since the fluctuation ((5r — 5s) 2 ) given in Eq.(26) becomes negative. The proba- 
bility distribution is not maximum there but represents a saddle point or a valley structure. On 



the other hand, new chiral fixed points bifurcate from r = s = \/ko/k2 when X± > 4ko/k 2 c 2 , 
as is described in Eq.(22). There are at most four chiral fixed points, depending on the rate 
coefficients. For k± = 0, the critical value of the quadratic rate coefficient k2 C is calculated as 

k 2c c 2 = 4fc (l + -^) 2 . (33) 
Ak 

When k 2 is larger than k2 C , the racemic fixed point is unstable, and correspondingly the final 
distribution Pf has double peaks associated with chiral states. 

When A — > 0, the final probability distribution is expected to vanish for Na > 0. On the 
line Na = 0, it is given by Pf(0,NR,N — Nr) in Eq.(17) but we are unsuccessful so far to 
obtain a closed analytic expression of the normalization factor N. When k\ = A = 0, Pf has 
peaks at chiral fixed points (r + , s_ = c — r + ) and (r_, s + = c — r_) with 



1 ± yi - 4fc /fc 2 c 2 

r± = c = ftp (34) 



The ee order parameter takes the value 



r± — s T 4k, 



±Jl--^L. (35) 

4. Numerical Integration of the Master Equation 

In this section we discuss on the time evolution of the probability distribution. It is ob- 
tained by integrating the master equation (10) by means of the fourth-order Runge-Kutta 
method. 31 Similarly to the previous section, three typical cases are discussed in the following 
subsections. For the numerical integration, we consider a system with the total number of 
active chemical species to be N = 100. 
4-1 Spontaneous production 

With a spontaneous production ko and a finite recycling as A = ko (and k\ = k 2 = 0), 
time evolution of the probability distribution is depicted in Fig. 1(a) by contour lines of a 
relative height interval with a quarter of the maximum height at several times. Whether 
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(c) (d) 

Fig. 1. Trajectories of probability contours obtained by numerically integrating the master equation 
(10) for spontaneous production (a) with (A = k ) and (c) without (A = 0) recycling, (b) and (d) 
represent final probability distributions with and without recycling, respectively. Other parameters 
are ki = k 2 = with N = 100. 

the system starts from a complete achiral state (Na, Nr, N$) = (100,0,0) or a chiral state 
(80, 20, 0), probability distributions keep a single-peak profile and converge to the unique final 
distribution. The final probability distribution shown in Fig. 1(b) agrees with the result given 
by Eq.(27). The peak position of the transient probability distribution follows the trajectory 
determined by the rate equations, drawn by continuous curves in Fig. 1(a). They approach a 
racemic fixed point N R = N$ = N/3 for k$ = A. 

Without recycling (A = 0), achiral substrate is consumed up ultimately, and the final state 
has a finite probability only at Na = or on a line Nr + N$ = N = 100. Also the remarkable 
fact is that without recycling the final state depends on the initial condition, as shown in 
Fig. 1(c); the probability distribution started from an initially achiral state Nr(0) = Ns(0) = 
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keeps its center on a racemic state Nr = N$ during the evolution, whereas that from a chiral 
initial state remains chiral even though the ee order parameter decreases. The final probability 
distribution given in the binomial form Eq.(30) is valid only for a racemic state, as shown in 
Fig. 1(d). 




Fig. 2. Trajectory of the probability contour for linear autocatalyses (ki = 10 _1 fco ) (a) with (A = fco 
) and (c) without (A = ) recycling, (b) and (d) represent final probability distributions with and 
without recycling, respectively. Other parameters are fca — with N — 100. 



4-2 Linear autocatalysis 

We now include a linear autocatalytic process with k\ = 0.1 fco- With a finite recycling 
as A = ko , the probability distribution converges to a unique profile with a single peak at 
a racemic state, as shown in Fig. 2(a): at the initial stage the system started from a chiral 
state with iV#(0) = 20 and Ng(0) = behaves differently from that started from an achiral 
state, but ultimately it turns to converge to a stable racemic fixed point. The final probability 
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distribution is the one given in Eq.(17), as shown in Fig. 2(b). 

When there is no recycling, the final probability depends on the initial condition, as shown 
in Fig. 2(c). The final probability staring from the achiral state Nr = N$ = is well described 
by the profile Eq.(32) determined by the detailed balance condition, as shown in Fig. 2(d). 




(e) (f) 

Fig. 3. Trajectory of the probability contour for nonlinear autocatalysis (k2 = 10~ 3 fco) (a), (c) with 
(A = ko) and (e) without (A = 0) recycling. Initial state for (a) is achiral, and for (c) is chiral. 
(b),(d) and (f) represent the final probability distributions corresponding to (a), (c) and (e), re- 
spectively. In (e) the final distribution obtained by numerical integration, represented by symbols, 
is quite different from the expected curve from the detailed balance condition, Eq.(17). Other 
parameters are k\ = with N = 100. 
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4-3 Quadratic autocatalysis 

With a quadratic autocatalysis, the probability distribution behaves qualitatively differ- 
ently from the previous cases. Even with a recycling process, the probability profile depends 
on the initial condition in the sense explained below. If the system starts from a completely 
achiral state without chiral species, N R = N$ = 0, the probability distribution is symmetric 
due to the symmetric dynamics, but the initial single peak splits around the unstable racemic 
fixed point, and the double-peak structure develops slowly, as shown in the contour evolution 
in Fig. 3(a). The parameters chosen are K2 = 10~ 3 fcoi K i = and A = k$. The final probabil- 
ity has a double-peak as shown in Fig. 3(b), and the profile is described by Eq.(17). On the 
contrary, when the system starts from a chiral state, as shown in Fig. 3(c), the probability 
distribution approaches a single peak profile, as shown in Fig. 3(d); the peak position corre- 
sponds to one of the double-peak structure in the symmetric case in Fig. 3(a). The probability 
distribution accompanies a long tail in the direction of the other peak position, but the sec- 
ond peak does not develop during the numerical integration upto the time of order 6/ A. This 
resembles to the symmetry breaking phase transition in equilibrium systems, and we may call 
that the system undergoes a chiral symmetry breaking. 

Without recycling (A = 0), the evolution of the system stops when the whole substrate 
molecules A are consumed N R + N$ = N, as shown in Fig. 3(e), quite similar to the previous 
two cases. However, the final probability started from an achiral state, shown by symbols 
in Fig. 3(f), is completely different from the A — > limit of Eq.(17), which is drawn by a 
continuous curve. It is not impossible, since at A = we cannot rely on the detailed bal- 
ance condition between the production and the recycling reaction. Then, what is the final 
probability distribution? We shall discuss this problem in the second next section. 

5. Eigenvalue Analysis of the Master Equation 

Before describing the probability distribution of a system without recycling, we consider 
the symmetry breaking for a nonlinear autocatalytic system with recycling in terms of the 
degeneracy in eigenvalues of the master equation. 

The master equation (10) is a linear differential equation for the probability distribution, 
rewritten by using an evolution matrix M as 



Here the evolution matrix M is related to the transition probabilities W, defined by Eqs.(7) 



dP(N A ,N R ,N s ;t) 



(N A ,N R ,N s \M\N' A ,N' R ,N' s )P(N' A ,N' R ,N' s -t). (36) 

N' A ,N' R ,N' S 



dt 



and (9), as 



(N A ,N R ,N S \M\N A ,N R ,N' S ) 
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W(N' A , N' R , N> s - N A , N R , N s ) for (N' A ,N' R , N' s ) ? (N A , N R , N s ), 

£ W(N A ,N R ,N S ^N A1 ,N R1 ,N S1 ) for (N' A , N' R , N' 3 ) = (N A ,N R ,N S ) (37) 

N A i,N m ,N sl 

and the matrix element is nonvanishing between states (N A , N R , N$) and (N A ,N R ,N' S ) with 
differences of numbers of molecules at most by one and N A + N R + Ng = N' A + N' R + N' s = N. 
By the use of eigenfunctions Vtj and eigenvalues A, of the matrix M as given by 

= Ai*i (38) 

the time development of the probability distribution is expressed by a series as 

oo 

P(N A , N R , N s ; t) = ^ a ie A ^. (39) 

i=0 

Since the probability distribution satisfies the conservation of the probability as ^2n r at s P = 
1, all the (A^+l)(A^+2)/2 eigenvalues Aj must be non-positive. Therefore, we order eigenvalues 
as Ao = > Ai > • • • Aj > Aj + i • • • . An equilibrium distribution Pf , if exists, is to be given 
by the eigenstate of the nondegenerate eigenvalue Ao = 0. We analyze three typical cases 
again. 

5.1 Spontaneous production 

In this case, the matrix M is simplified by introducing matrices A, A*, R, R*, S, S* 
which have mostly zero elements except the following nonzero ones as 

(N A - 1,N R ,N S \A\N A ,N R ,N S ) = N A , (N A + 1, N R , N S \A*\N A , N R , N s ) = 1, 
(N A , N R - 1, iV s |R|AU, N R , N s ) = N R , (N A , N R + 1, N S \R*\N A , N R , N s ) = 1, 
(N A ,N R ,Ns-l\S\N A ,N R ,N s ) = N s , (N A , N R , N s + 1\S*\N A , N R , N s ) = 1. (40) 
These matrices satisfy the following nonzero commutation relations 

[A, A*] = [R, R*] = [S, S*] = 1 (41) 

and other commutators vanish. Regarding the state (0,0,0) as a vacuum, we obtain 
(A*)^ |0, 0, 0) = \N A , 0, 0), and thus A* may be regarded as an normalized creation operator 
of A molecules. Similar interpretation holds for other matrices, and thus three bilinear prod- 
ucts A*A, R*R, S*S represent number operators. The time-evolution matrix M is expressed 
as a bilinear form of these matrices as 

M = (R* - A*)(A; A - AR) + (S* - A*)(k A - AS). (42) 
By introducing new combinations defined as 

C = A + R + S. cj = *V + Wl* + W 



Ci — R — S, 



2k + A 
R* - S* 

2 ' 



OA* I T> * I 

C 2 = -2A: A + AR + AS, C* = -±- -f , (43) 

4/eg + zA 
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one can prove that these new matrices satisfy the commutation relations similar to Eq.(41) 
as [Cj, C£] = Sij, [Cj, Cj] = [C*,C*] = 0, and the time evolution matrix is represented in a 
diagonal form as 

M = -AC*Ci - {2k + A)C*C 2 . (44) 

M's independence of CgCo is related to the total number conservation in the present system. 
Since the eigenvalues of CJCi and C2C2 are non-negative integers n\, ri2 = 0, 1, 2 ■ ■ ■ , and 
n i + n 2 < A/", the eigenvalues of the time evolution is written as 

A ni;n2 = -Am - (2k + A)n 2 . (45) 

The zero eigenvalue Aoo is nondegenerate as long as A > 0, and the system relaxes to a unique 
final state: No phase transition is expected. 

On the other hand, without recycling A = there is a large degeneracy in the eigenvalues. 
Especially for a zero eigenvalue there are N + 1 degeneracy since A nij o = for any n\ from 
to N. It represents the non-ergodicity of the system and the final state depends on the initial 
condition, as is anticipated from the rate equation analysis. 

5.2 Linear autocatalysis 

With a linear autocatalysis, we have failed so far to obtain analytic expressions of eigen- 
values of the time-evolution matrix M, and we have to rely on the numerical method. We 
use the subroutine "dgeev" in LAPACK to obtain eigenvalues and eigenfunctions of a gen- 
eral matrix, but because of the algorithmic restriction, the maximum size of N is limited to 
50. Because of the smaller size N in this section than the one in the previous sections, the 
parameters are chosen to be a little larger as K\ = 0.2/co, ki = 0, A = ko- First few largest 
eigenvalues are depicted in Fig. 4(a) as a function of the total number of chemical reactants N. 
With a finite recycling A > 0, the second largest eigenvalue Ai increases with N, but seems 
to saturate before it approaches Ao = 0. This finite gap between Ai and means the absence 
of non-ergodicity and of the phase transition. The 0th eigenfunction Fig. 4(b) looks similar to 
the final probability distribution for a larger system with N = 100, shown in Fig. 2(b). The 
1st eigenfunction in Fig. 4(c) is assymmetric around Nr = N$ line. 

5.3 Quadratic autocatalysis 

In the case with a quadratic autocatalysis with parameters, K2 = 4 x 10 _3 A:o and A = 
ko, k\ = 0, the second largest eigenvalue Ai approaches within a system size upto A?" =50, 
whereas the third largest A2 saturates to a finite value, as shown in Fig. 5(a). According to 
Eq.(33), the critical size for the chiral symmetry breaking is calculated to be N c = 2(1 + 
A/ 'ko)(n2/ ko) 1 ^ 2 ~ 39.5 in the present choice of the parameters. In a semi-logarithmic plot of 
the eigenvalues in Fig. 5(b), the splitting of — Ai and — A2 is in fact taking place close to this 
critical size N c . Near the maximum possible size A?" = 50, the size-depence of the eigenvalue 
Ai looks like to be exponential as Ai = —2 x 10 3 exp(— 0.23N). Of course, the eigenvalue 
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Fig. 4. (a) Eigenvalues of the time-evolution matrix with a spontaneous and a linear autocatalytic 
production and recycling; m = 0.2fco, A = fco and ki = 0. (b) The cigcnfunction of the zero 
eigenvalue Ao = 0, and (c) <I>i of the second largest eigenvalue Ai at N — 50. 

degeneracy seems just have commenced, and the true asymptotics might be accomplished 
for much larger systems. Therefore, the formula is only tentative. Still, the emergence of 
non-ergodicity exists for sure, and the phase transition takes place in a thermodynamic limit 
N ^ co. 

We know from the analysis in §2 that the final distribution with this autocatalysis has a 
symmetric double peak structure. That corresponds to the eigenfunction Ior the eigenvalue 
Ao = 0, shown in Fig. 5(c), and looks similar to the final probability distribution in Fig. 3(b) 
for a larger system N = 100. The eigenfunction corresponding to Ai is asymmetric, as 
shown in Fig. 5(d). If one starts from an initial state with a small chirality, as in the case 
of Fig. 3(c), some component of \&i is mixed in the initial state so as to enhance one peak 
and suppress the other peak in = Pf- As time passes, i > 2-components die off, but 
the two components (^/q? ^i) remain over the duration 1/Ai, and the asymptotic probability 
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distribution is single peaked at a chiral state, as shown in Fig. 3(d). 

When the recycling vanishes (A — ► 0), not only one but many eigenvalues approach zero, 
irrespective of the system size, and the system is completely non-ergodic. The final state 
depends on the initial state. We cannot apply the analysis employed here, and have to consider 
the problem in a completely different way in the next section. 
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Fig. 5. (a) A few largest eigenvalues of the time-evolution matrix with a spontaneous and a quadratic 
autocatalytic production and recycling; K2 = 4 x 10~ 3 fco, A = ko and fei = 0. The second largest 
eigenvalue Ai approaches as the total number N increases, (b) Semi-logarithmic plot of the 
magnitude of a few largest eigenvalues. The asymptotic behavior Ai w — 2 x 10 3 exp(— 0.23A) fits 
well the second largest eigenvalue Ai. (c) The eigenfunction with the zero eigenvalue Ao = 0, 
and (d) of the second largest eigenvalue Ax at N = 50. 



6. Directed Random Walk Model of Non- Recycling System 

So far there remains an apparent discrepancy between the final probability distribution 
for the nonrecycling system obtained by the time integration of the master equation and 
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that obtained as a nonrecycling limit of the analytic result derived by the detailed balance 
consideration. For a spontaneous and for a linear autocatalytic cases, these two analyses give 
the same results, whereas for a nonlinear autocatalytic case they are different. In fact, in the 
nonrecycling case, the detailed balance condition cannot be imposed from the beginning, since 
there is no back reaction to balance the production process. Therefore, one has to consider in 
a completely different way. 

Here we propose a toy model which may be relevant to the stochastic evolution of an 
autocatalytic system without recycling. It is a random walk model on a square lattice (Nr, Ng) 
in a triangular region < Nr, Ns, (Nr + N$) < N. A random walker can only make directed 
walks to the right or upwards: Nr —* Nr + 1 or N$ — > Ng + 1. The transition probability 
is Eq.(7) to the right, and the corresponding one to upwards. This means that a walker on a 
site (Nr, Ns) stays on the site for a waiting time interval 

r(N R , N S ) = [2kQ - ^ {Nr - - K2(Ar 2 + n2)]Na ( 46 ) 

with Na = N — Nr — Ns, and then jumps to the right and upwards with probabilities p r and 
p u given as 

k + Kl N R + k 2 N 2 r 



Pr(N R ,N s ) 



2k + Ki(N R + N S ) + k 2 (N 2 r + N 2 S ) ' 



Pu(Nr, N S ) = 2k) + Ki{Nr + ^ + + , (47) 

respectively. The average position of the walker is expected to follow the evolution similar to 
the rate equations Eq.(5) with A = 0. The remarkable point of this model is that though the 
waiting time t(Nr, Ns) gets longer as Na decreases, the rates of the transitions to the right 
and upwards do not depend on Na nor vanish when Na = 0. For a walker to reach on a line 
Nr + Ns = N, it may take an exceedingly long time for an external observer, but it requires 
only N steps for a walker himself. 

Let the walker start from the origin (Nr,Ns) = (0,0). After a time t^ 1 ' = r(0, 0) = 
l/2k()N he jumps to the right or upwards with the same probability, and lands on a line 
Nr + Ns = 1. Wherever he is on a line Nr + = 1, he makes a second jump at a time 
f(2) = f(!) + [(2& + Ki + k 2 )(N - l)]" 1 to the line N R + N s = 2. For a nonautocatalytic case 
with k 2 > 0, the time of the third ( and the further) jump depends on a position Nr and Ns, 
but when k 2 = the walker jumps to the next line Nr + Ns = 3 at a same time t^. In fact, 
after the n-th jump, the walker corresponding to the linear autocatalysis is somewhere on a 
line Nr + Ns = n, and jumps to the next line Nr + A^s = n + 1 at the same time independent 
of the position on the line. In this case, we can calculate the probability evolution analytically. 
6.1 Directed Random Walker corresponding to Spontaneous and Linear Autocatalysis 

A walker is making a random walk corresponding to the spontaneous production and 
a linear autocatalysis with a finite k$ and k\ but without nonlinear autocatalysis (k 2 = 0). 
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When he is on a lattice site (Nr, Ns), he has to wait for a next jump a time interval t(Nr, N$) 
which depends only on the sum Nr + Ns = n but not on Nr and Ns individually, so that 
the waiting time can be denoted by r n . This sum n also represents the number of jumps 
he has made, and totally t^ = Y17=o r * time has elapsed since he started from the origin 
(Nr, Ns) = (0, 0) till he reaches the present site. After n jumps, the walker is somewhere on a 
line Nr + Ns = n, and we denote the probability as Prw(Nr, Ns;n). It is non-zero on a line 
Nr + Ns = n but zero anywhere else. He should have come from the left or from the down. 
Therefore, the probability should satisfy the relation 

P RW {N R ,Ns;n) 

= Pr {N R - 1, N s )Prw{Nr - 1, N s ;n - 1) + Pu (N R , N s - 1)Prw{Nr, N s -l;n-l) 

= (kg + m(N R - 1))Pr W (Nr - 1, N s ; n - 1) + (kg + Kl (N s - 1))Prw(Nr, N s -l;n-l) 

(2k + Kin) 

+ N S )i ngr^o + m n?sH*o + ^) p rnnn , UR , 

= ^!Ay Ul R +Ns ~\2k + Kl l) ^ 

The combination factor represents the number of ways of arranging the order of right and 
upwards jumps among the Nr + Ns = n total jumps from the initial state at the origin: 
Prw(0, 0; 0) = 1. After n = N jumps, the walker reaches to the final line Nr + Ns = N, and 
the distribution of the walker Prw(Nr, Ns) agrees with the final distribution (32) discussed 
in §3.2 in the limit of A — ► 0, and naturally with the result of time integration, shown in Fig. 
1(d) and 2(d). 

If the walker has started from an arbitrary site (Nrq, Nso), he can be on a site Nr > Nrq 
and Ns > Nso with a probability 

PmviNR - Ns)= (N R -N m y.(N s -N S0 y. iSg^ ' (49> 

It explains the asymmetric final probability distribution obtained by starting from a chiral 
initial state, shown in Fig. 2(c). 

6.2 Directed Random Walker corresponding to Nonlinear Autocatalysis 

The same analysis cannot be performed on the time evolution of the random walker cor- 
responding to the nonlinear autocatalysis, since his dwelling time on a site (Nr, Ns) depends 
not only on the sum Nr + Ns but also Nr and Ns individually. But after the nth step, the 
walker is on a line A^^ + Ns = n. Only the arrival time depends on the path he took from 
the start at (0,0). If he comes on the site after the nth jump, he has come from the left 
(Nr — l,Ns) or from below (Nr,Ns — 1). Therefore, the total probability Prw(Nr, Ns : n) 
that he passes the site Nr + Ns = n satisfies the relation Eq.(48); 

Prw(Nr, N s ; n) = Pr (N R - 1, N s )Pr W (Nr - 1, N s ; n-1) 

+ Pu (Nr,N s - 1)Pr W (Nr,N s - l;n - 1). (50) 
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Because the denominators of p t (Nr — 1, N$) and p u (Nr, N$ — 1) are not common, we cannot 
calculate the distribution analytically. However, by starting from Prw(Nr = 0,Ns = : 
0) = 1, we can calculate Prw(Nr, N$ : n) for any combinations of Nr and Ng. Since the 
random walker passes the particular site (Nr, N$) once and for ever, and he stops on a line 
N = Nr_ + Ns if one considers the time evolution, the final probability distribution of the 
walker is given by Prw(Nr, N$; N). 

For the nonlinear autocatalytic system of a size TV = 100 with reaction parameters K2 = 
10 _3 fco and k\ = A = 0, we obtain the final probability distribution which has a single 
maximum at a racemic state, in agreement with the result of time integration, as shown in 
Fig. 6(a). 




N R K 2 /k) 



(d) (e) 

Fig. 6. (a) Final probability distribution of a nonlinear autocatalytic system: N = 100, = 
10 _3 fco, Ki = A = 0. Symbols are the result by the time integration of the master equation and a 
continuous curve represents numerical calculation of the directed random walk model, (b) A trace 
of the probability maxima for K2 = 10~ 3 . (c) Variation of the final probability distribution as a 
function of N. (d) A trace of the probability maxima for «2 = 10~ 2 . (e) Critical number N c for 
symmetry breaking at various strength of quadratic autocatalytic rate constant K2- The fit line is 
N c = 1.33( K2 /fc )-°- 75 . 

In Fig. 6(a) the probability distribution has a single peak, but in fact, the profile depends 
on the total number of chemical reactants N. The peak positions of the probability distribution 
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Prw are shown in Fig. 6(b) up to a size N = 400. It has a single peak for small Nr + N$ = 
N < N c , whereas for large sizes N > N c it has symmetric double peaks. The probability start 
to split to have a double-peak structure around N c ~ 235, and for N=400 the probability has 
clear double peaks, as shown in Fig. 6(c). The critical value in terms of the concentration is 
&2 C c 2 = H2N 2 ~ 55.2&o, quite larger than the value estimated from the A — > limit value 
&2 C c 2 = 4&o- Another fact to be noticed in Fig. 6(b) is that the traces of the two peaks 
are asymptotically parallel to two axes: Nr at the peak position, for example, is fixed, but 
Ns = N — Nr increases with N, and so does the magnitude of the ee order parameter 
\</>\ = \N R - N s \/N. 

The critical number N c for double peaks actually depends on ^2/^0- F° r example, on 
increasing K2 to K2 = 10 _2 &;o, the traces of probability maxima is scaled down, as shown in 
Fig. 6(d). The critical size for the double-peak structure in this case is N c ~ 43, and in terms 
of the rate constant k<i its critical value is A;2 C c 2 = K2N 2 = 18.5/co- By varying ^2/^0 from 
10 -4 to 10 _1 , one finds that the critical size N c depends on K2/&0 as N c = 1.33(k2//co) -0 ' 75 , 
as shown in Fig. 6(e). We will show later that the exponent —3/4 has a certain meaning. 

When K2 is close to k$ and N c is of the order unity, the discreteness of the integer number 
comes into play in the formula for iV c . In fact, at K2 ~ ^0 the probability peak takes place 
on the boundary, Nr = and N$ = 0. A similar behavior is already found for the linearly 
autocatalytic case. This represents a single mother scenario of the homochirality by the au- 
tocatalytic reactions: When one starts from a completely achiral state without any chiral 
ingredients initially, the first chiral species produced randomly by the spontaneous reaction 
converts all the achiral substrate A into her descendants before the second mother of different 
chirality is born. 

We now study the case when K2 is much smaller than fcoj K 2 *C &o- Even in this case, the 
probability distribution has double peaks when the total number of reactants is large enough, 
as shown in Figs.6(b)-(d). By starting from an achiral state without any chiral species, the ini- 
tial production is governed by the spontaneous production. Only after n steps of spontaneous 
production where 

K 2 n 2 ~ k , (51) 

there are enough chiral products for quadratic autocatalysis: n ~ \/fco/ K 2 ^ 1- Here the 
probability distribution is Gaussian centered at a racemic state Nrq = N$o = re/2 with a 
width n. In terms of the difference xq = Nrq — N$o at a step n = Nrq + AgO) the probability 
is approximately given as 

P(x )dx oc e~ x ° /2n dx . (52) 

After the initial stage of n steps, the quadratic autocatalysis comes into play. The state 
with (Nrq, Nso) evolves deterministically according to the rate equations of the quadratic 
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autocatalysis. At the N-th step a state reaches to (Nr, N$) with Nr + Ns = N and related 
to the "initial" state (Nrq, Nso) as 

J— L = J— L (53) 

Nr N s N m jV s „ 1 ' 

or in terms of the chirality x = Nr — Ns as 

n 2 -xl n 2 / N 2 -n 2 2 \ 
X0 = X N^*Nt X { 1 + ^^ X+ ---)- (54) 

The second approximation is valid close to the racemic state with |xo/n|, \x/N\ <C 1. The 

probability distribution of x around a racemic state \x/N\ <C 1 is obtained as 

P(x)dx = P(x Q )dx oc e -* 2 °/ 2n ^dx ~ e- Ax2 ' 2+ °^dx (55) 



with the coefficient A written as 



^=^"4 (56) 



AT4 TV 2 

for N S> n S> 1. The coefficient ^4 becomes negative for 

indicating that the probability distribution is minimum at x = 0. It means that the probability 
distribution has peaks at somewhere different from the racemic point x = 0. It should be a 
double peak structure because of the symmetry. The power law dependence Eq.(57) of the 
minimal number of chemical reactants N on the ratio of rate constants with an exponent -3/4 
is in good agreement of the numerical data, shown in Fig. 6(e). 

This mathematical analysis is qualitatively interpreted in the following manner. With a 
small K2, the chiral species are produced spontaneously in the initial stage and the probability 
distribution has a racemic single peak. After the numbers of chiral species n are large enough as 



n ~ \/ko/K2, the quadratic autocatalysis increases the population of the majority enantiomer 
more than that of the minority ones, and when there are totally chiral molecules, the 
probability of the racemic state decreases with a fluctuation about A 4 /n 3 larger than the 
Gaussian width N. Furthermore, the nonlinear relation between the initial enantiomeric excess 
or xo to the final one x changes the Jacobian. Therefore, after a sufficient steps of nonlinear 
growth N ~ (fco/ K 2)~ 3 ^ 4 ) probability distribution takes a double peak structure. The rate 
coefficient for the quadratic autocatalysis should be large as K2A 4 / 3 > &o f° r the probability 
distribution to have a double-peak structure, but this condition is weaker than that for the 
single mother scenario; K2 > ko- 

7. Summary and Discussions 

The probability distribution of the populations of chiral species in the reaction of chiral 
molecule production is studied by means of stochastic master equation. With a recycling back 
reaction, the system relaxes to a unique final distribution, which can be obtained analytically 
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by assuming a detailed balance condition. The final probability distribution has a racemic sin- 
gle peak for systems with a spontaneous and a linear autocatalytic reactions. With a quadratic 
autocatalysis, it has a double-peak structure, indicating the chiral symmetry breaking. 

Without recycling, the probability is shown to depend on the initial condition. By starting 
from a complete achiral initial state without chiral ingredients, the final distribution agrees 
with that predicted by the weak recycling limit obtained by the detailed balance condition 
in a spontaneous and a linear autocatalytic cases, but in a quadratic autocatalytic system 
detailed balance condition leads to erroneous final probability different from that obtained 
by the time-integration of the master equation. Without recycling, in fact, there is no reason 
that the detailed balance condition holds. Instead, we presented another directed random 
walk model for the non-recycling system. The final probability is explained in all the cases 
without recycling by this model. With the quadratic autocatalysis, the probability develops 
double peaks as the rate of quadratic autocatalysis k<i increases. Quadratic autocatalysis not 
only broadens the fluctuation produced initially by spontaneous production but also increases 
the density of states away from the racemic state. Therefore, if the rate coefficient ki is larger 
than that of the spontaneous production ko by a factor about V 2 /N 4: ^ 3 , double peaks develop 
in the probability distribution, i.e. if the reaction takes place in a small volume V with a large 
number of reactants N. 

In the limit of a very small rate of the spontaneous production ko <C k\/V or ko <C k^jV 2 ^ 
the probability distribution has double peaks where one of the enantiomer is missing. This 
corresponds to the single mother, or chirality Eve, scenario for the homochirality. After the first 
chiral molecule is produced spontaneously and randomly, all the available achiral substrate 
molecules are converted to the first mother's type by linear or quadratic autocatalysis, before 
the second mother of another enantiomeric type is born spontaneously. This rare situation 
cannot be described by the rate equation but only by the stochastic method. 

As for the double peaks in the probability distribution found in experiments of the Soai 
reaction, they are not sharp ones at the perfect homochirality. Therefore, the single mother 
scenario is improbable, and the quadratic autocatalysis with a moderate strength seems most 
plausible. However, if the single mother scenario is accompanied with some imperfections as 
epimerization process 32 or erroneous catalysis, peaks may shift to smaller ee values. More 
careful study is necessary. 
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